function [Mint] = distribute_proportionally_2_integers(Mtotals, weights)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Suppose I have Mtotal, which is an integer number. I need to distribute it proportionally
	% into K bins, taking into account the weights r_k. But the results need to be integers.
	%
	% This function does just that, respecting the total at the end.
	% It first computes the theoretical (non-integer) shares. Then it rounds off the values into integers
	% Finally, it figures out whether the total is less or more than the original total, and does the rebalancing
	% by identifying the bins that have been the most favored (or disadvantaged) by the rounding.
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% Mtotal:			1 x dim2  (integers)
	% weights:			K x 1 or K x dim2
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% Mint:				K x dim2 (integers such that sum(M) = Mtotal )
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% 1. Compute the theoretical (non-integer) shares
	M = (weights./sum(weights)) .* Mtotals; % K x dim2
	
	% 2) Round of the values into integers
	Mint = round(M);

	% 3) Compute the rounding errors and the total errors
	round_errors = Mint - M;
	diffs = sum(Mint) - Mtotals; % 1 x dim2
	
	% 4) Fix the errors
	for jj = 1:size(Mint,2)
		diff_jj = diffs(jj);
		if diff_jj < 0 % I still have a positive number diff_jj that I need to distribute into the K bins.
			% Find the smallest rounding errors.
			[~, idxes] = mink(round_errors(:,jj),abs(diff_jj));
			Mint(idxes,jj) = Mint(idxes,jj) + 1;
		end
		if diff_jj > 0 % I have a positive number diff_jj that I need to remove from the K bins.
			% Find the biggest rounding errors.
			[~, idxes] = mink(-round_errors(:,jj),abs(diff_jj));
			Mint(idxes,jj) = Mint(idxes,jj) - 1;
		end
	end
end
